Serie de Tiempo

¿Que son?

Una serie de tiempos, a diferencia de cualquier otro data set, tiene como característica un campo con información temporal u orden. Muchas disciplinas, como finanzas, administración pública, energía, retail, salud, están dominadas por información de serie de tiempos.

Demos algunos ejemplos: 1) Cotizaciones de acciones. 2) El INDEC calculará el producto interno bruto del país sobre una base anual. 3) 2020, año COVID seguro han seguido y visto los casos semanales de infecciones. 4) Otro ejemplo es el calentamiento global en la era post industrial.

Podemos decir que una serie temporal o cronológica es una sucesión de datos medidos en determinados momentos y ordenados cronológicamente

Algo importante en series temporales es la autocorrelación: Los valores de una serie temporal está correlacionada con valores anteriores. Si no existe ningún tipo de correlación, entonces fenómeno no puede ser modelado como una serie de tiempo.

Hablemos ahora de sus características:

  • Tendencia secular o regular, indica la marcha general y persistente del fenómeno observado, es una componente de la serie que refleja la evolución a largo plazo. Por ejemplo, el uso creciente de Internet en la sociedad, independientemente de que en un mes concreto en un país, por determinadas causas, haya una baja en la utilización de Internet.

  • Variación estacional o variación cíclica regular, el movimiento periódico de corto plazo. Se trata de una componente causal debida a la influencia de ciertos fenómenos que se repiten de manera periódica en un año (las estaciones), una semana (los fines de semana) o un día (las horas puntas) o cualquier otro periodo. Recoge las oscilaciones que se producen en esos períodos de repetición. Por ejemplo, el tráfico en la autopista que reporta picos durante determinadas franjas horarias y días

  • Variación cíclica, el componente de la serie que recoge las oscilaciones periódicas de amplitud superior a un año. Suelen deberse a la alternancia de etapas de prosperidad económica (crestas) con etapas de depresión (valles).Por ejemplo, ciclos económicos, recesiones

  • Variación aleatoria o ruido, accidental, de carácter errático, también denominada residuo, no muestran ninguna regularidad (salvo las regularidades estadísticas), debidos a fenómenos de carácter ocasional. Por ejemplo tormentas, terremotos, inundaciones, huelgas, guerras, avances tecnológicos, etcétera.

Usaremos en esta notebook la libreria astsa y sus datasets.

library(astsa) #cargo libreria

Aquí tenemos la serie de ganancias trimestrales por acción de Johnson & Johnson. Tiene algunas características comunes de los datos de series de tiempo, tendencia al alza, estacionalidad en el sentido de que el segundo y tercer trimestres generalmente hacia arriba, mientras que el cuarto trimestre generalmente es hacia abajo. Además, existe heterocedasticidad porque, a medida que aumenta el valor del activo, los pequeños cambios porcentuales se convierten en grandes cambios absolutos.

#cargamos la base jj
data(jj)
plot(jj, main = "Johnson & Johnson Quarterly Earnings per Share", type = "c")
text(jj, labels = 1:4, col = 1:4)

La segunda serie son las desviaciones anuales de la temperatura global. Los datos son desviaciones de la temperatura promedio entre 1960 y 1980. Los datos tienen una tendencia generalmente positiva, pero la tendencia no siempre es positiva. A diferencia de los datos de Johnson y Johnson, esta serie no tiene un componente estacional y es homoscedástica.

data("globtemp")
plot(globtemp, main = "Global Temperature Deviations", type= "o")

La tercera serie son los rendimientos semanales del S&P 500 (índice bursátil estadounidense basado en 500 grandes corporaciones). Las devoluciones son el cambio porcentual por período de tiempo. A diferencia de las otras series, esta serie no tiene tendencia ni estacionalidad. De hecho, parece que no hay ningún patrón en la serie (excepto que de vez en cuando, la varianza es grande). Este es un ejemplo de un tipo particular de proceso llamado ruido.

data(xts)
## Warning in data(xts): data set 'xts' not found
data("sp500w")
plot(sp500w, main = "S&P 500 Weekly Returns")

Descomponiendo series de tiempo

Una serie de tiempo la podemos descomponer: tendencia, ciclo y estacionalidad. Sacando esto del modelo, solo queda el residuo el cual, en un buen modelo, no debe tener correlación serial ya que el modelo debería haber capturado el patrón. A este ruido blanco gaussiano sin correlación serial se lo denomina comúnmente como Random Walk, siendo el valor previo de la serie más ruido blanco que es aleatorio.

descomposejj <- decompose(jj, type = c("additive", "multiplicative"), filter = NULL) #uso de descompose en r para descomponer jj 
plot(descomposejj) #dibujamos la descomposicion

Para trabajar con series de tiempo necesitamos remover la tendencia y forzar los datos a ser estacionales a esto se lo denomina DIFERENCIACION. La diferenciación analiza la diferencia entre el valor de una serie de tiempo en un momento determinado y su valor anterior.

Veamos a continuacion una tendencia y su diferenciaciòn.

plot(jj)
title("JJ trend stationary")

jj_diff <- diff(jj)
plot(jj_diff)
title("JJ random walk")

# ARMA

Cualquier serie de tiempo estacionaria se puede escribir como una combinación lineal de ruido blanco. El modelo ARMA tiene esta forma, por lo que es una buena opción para modelar series de tiempo estacionarias.

Los modelos ARMA son modelos de regresión de series de tiempo. Si recordamos, en la regresión tenemos una variable dependiente (Y), una variable independiente (X) y regresa linealmente Y sobre X. Una suposición crucial es que los errores son independientes, normales y homocedásticos. En otras palabras, los errores son ruido blanco. [WN] El ruido blanco es una secuencia de normales independientes con varianza común. Eventualmente verá que los modelos de series de tiempo se construyen alrededor del ruido blanco.

[AR]

[AR] Con las series de tiempo, puede hacer una regresión de hoy a ayer, y esto se denomina regresión automática (o autoregresión). O sea predecimos sobre el mismo valor del ciclo anterior. En este caso, lo que sucede hoy es la variable dependiente y lo que sucedió ayer es la variable independiente.

[MA]

[MA] Por lo general, los datos de las series de tiempo están correlacionados, y asumir que los errores no están correlacionados puede conducir a malos pronósticos. Una forma de superar el problema es utilizar una media móvil para los errores. MA realiza cálculos basados en el ruido en los datos junto con la pendiente de los datos Es también un suavizador de los valores en la medida en que el n se hace más grande. Esta media móvil también puede ser ponderada según el período para darle mayor fuerza a los valores más cercanos al tiempo presente.

Generacion de modelos ARMA

Vamos a ir generando datos a partir de varios modelos ARMA. Generando 200 observaciones y graficando el resultado con el modelos arima.sim de r.

# Generamos ruido blanco
WN <- arima.sim(model=list(order = c(0,0,0)),n=200)
plot(WN)

# Generamos modelo AR(2) con parametros 1.5 and -.75
AR <- arima.sim(model=list(order = c(2,0,0), ar=c(1.5,-0.75)),n=200)
plot(AR)

# Generamos modelo MA(1) con parametro .9 
MA <- arima.sim(model=list(order = c(0,0,1), ma=0.9),n=200)
plot(MA)

Viendo estos graficos vemos que no es tan sencillo identificar modelos ARMA. Vimos el modelo AR y el modelo MA y resultan bastante similares, no pudiendo identificarlos simplemente mirando los datos.

Identificaciòn ARMA (ACF y PACF)

Las herramientas que se utilizan para identificar los órdenes del modelo son la función de autocorrelación (o ACF) y la función de autocorrelación parcial (o PACF). Si un proceso es AR puro, entonces el ACF se reducirá y el PACF se interrumpirá en el lag p. Para un MA puro, es lo opuesto: el PACF se reduce y el ACF se corta en el lag q. Si ambos están disminuyendo, entonces el modelo es ARMA.

Utilizaremos sarima() del paquete astsa para ajustar fácilmente los modelos a los datos. El comando produce un gráfico de diagnóstico residual que se puede ignorar por el momento (lo veremos mas adelante).

ACF - PACF - [AR]

Generamos 100 observaciones en un modelo AR(1)

# Generamos 100 observaciones en un modelo AR(1) 
ar <- arima.sim(model = list(order = c(1, 0, 0), ar = .9), n = 100) 

# Plot 
plot(ar)

Evaluamos ahora el ACF y PACF. Podemos observar el ACF y como se reduce y el PACF que se interrumpe en el lag p=1.

# Dibujamos ACF y PACF 
plot(acf2(ar))

# Fit del modelo AR(1) y examinamos la t-table

#ar_fit<-sarima(ar, p=1,d=0,q=0)
#ar_fit$ttable

ACF - PACF - [MA]

Hagamos lo mismo pero ahora con un modelo MA.

# Generamos 100 observaciones en un modelo MA(1) 
ma <- arima.sim(model = list(order = c(0, 0, 1), ma = -.8), n = 100)
plot(ma) #plot

Para un MA puro, el PACF se reduce y el ACF se corta en el lag q =1.

# Dibujamos ACF y PACF 
plot(acf2(ma))

# Fit modelo MA(1) y t-table
#ma_fit<-sarima(ma, 0,0,1)
#ma_fit$ttable

ACF - PACF - [ARMA]

Generamos ahora modelos ARMA(2,1)

# Generamos 250 observaciones en un modelo ARMA(2,0,1) 
x_arma <- arima.sim(model = list(order = c(2, 0, 1), ar = c(1, -.9), ma = .8), n = 250)
# Plot 
plot(x_arma)

Observemos ahora el PACF y ACF para un modelo ARMA. Vemos que ambos están disminuyendo, entonces el modelo es ARMA.

# Plot PACF - CFA
plot(acf2(x_arma))

# Fit modelo ARMA(2,0,1) y t-table
#x_fit<-sarima(x,2,0,1)
#x_fit$ttable

Evaluaciòn de un modelo

Para analizar el mejor modelos usaremos dos métricas:El AIC y el BIC así como el análisis de residuos.

AIC y BIC

El criterio de información de Akaike (AIC) es una medida de la calidad relativa de un modelo estadístico, para un conjunto dado de datos. Como tal, el AIC proporciona un medio para la selección del modelo. En estadística, el criterio de información bayesiano (BIC) o el más general criterio de Schwarz (SBC también, SBIC) es un criterio para la selección de modelos entre un conjunto finito de modelos.

El AIC tiene un k=2 y el BIC tien un k=log(n)

El objetivo es encontrar el menor BIC o AIC.

Analisis de los residuos

Si analizamos gráficamente los modelos debemos tener en cuenta

  • Los residuos estandarizados deben comportarse como una secuencia de ruido blanco con media cero y una varianza. Examinanmos los gráficos de residuos en busca de desviaciones de este comportamiento.

*El ACF de muestra de los residuos debe verse como el del ruido blanco. Examinamos el ACF en busca de desviaciones de este comportamiento.

*La normalidad es un supuesto esencial al instalar modelos ARMA. Examinamos la gráfica Q-Q para ver si hay desviaciones de la normalidad e identificar valores atípicos.

*Utilizamos la gráfica de estadística Q para ayudar a probar las desviaciones de la blancura de los residuos.

BIC & AIC + Analisis de residuos - Ejemplo

Vamos a trabajar sobre la base de datos “varve” de r comenzando con un modelo MA1 e incrementando el orden para entender si el modelo mejora o no.

Esta base “Varve” contiene depósitos sedimentarios de un lugar en Massachusetts durante 634 años, comenzando hace casi 12.000 años.

data("varve") #cargo la base de varve
?varve #Depósitos sedimentarios de un lugar en Massachusetts durante 634 años, comenzando hace casi 12.000 años.
dl_varve <- diff(log(varve)) #aplicamos logaritmo a la base
# Fit MA(1) a dl_varve.   
MA_fit <- sarima(dl_varve,0,0,1)
## initial  value -0.551780 
## iter   2 value -0.671633
## iter   3 value -0.706234
## iter   4 value -0.707586
## iter   5 value -0.718543
## iter   6 value -0.719692
## iter   7 value -0.721967
## iter   8 value -0.722970
## iter   9 value -0.723231
## iter  10 value -0.723247
## iter  11 value -0.723248
## iter  12 value -0.723248
## iter  12 value -0.723248
## iter  12 value -0.723248
## final  value -0.723248 
## converged
## initial  value -0.722762 
## iter   2 value -0.722764
## iter   3 value -0.722764
## iter   4 value -0.722765
## iter   4 value -0.722765
## iter   4 value -0.722765
## final  value -0.722765 
## converged

#MA_fit$ttable

# # Fit MA(2) a dl_varve. 
MA2_fit <- sarima(dl_varve,0,0,2)
## initial  value -0.551780 
## iter   2 value -0.679736
## iter   3 value -0.728605
## iter   4 value -0.734640
## iter   5 value -0.735449
## iter   6 value -0.735979
## iter   7 value -0.736015
## iter   8 value -0.736059
## iter   9 value -0.736060
## iter  10 value -0.736060
## iter  11 value -0.736061
## iter  12 value -0.736061
## iter  12 value -0.736061
## iter  12 value -0.736061
## final  value -0.736061 
## converged
## initial  value -0.735372 
## iter   2 value -0.735378
## iter   3 value -0.735379
## iter   4 value -0.735379
## iter   4 value -0.735379
## iter   4 value -0.735379
## final  value -0.735379 
## converged

#MA2_fit$ttable

# Fit ARMA(1,1) a dl_varve
ARMA_fit <- sarima(dl_varve,1,0,1)
## initial  value -0.550994 
## iter   2 value -0.648962
## iter   3 value -0.676965
## iter   4 value -0.699167
## iter   5 value -0.724554
## iter   6 value -0.726719
## iter   7 value -0.729066
## iter   8 value -0.731976
## iter   9 value -0.734235
## iter  10 value -0.735969
## iter  11 value -0.736410
## iter  12 value -0.737045
## iter  13 value -0.737600
## iter  14 value -0.737641
## iter  15 value -0.737643
## iter  16 value -0.737643
## iter  17 value -0.737643
## iter  18 value -0.737643
## iter  18 value -0.737643
## iter  18 value -0.737643
## final  value -0.737643 
## converged
## initial  value -0.737522 
## iter   2 value -0.737527
## iter   3 value -0.737528
## iter   4 value -0.737529
## iter   5 value -0.737530
## iter   5 value -0.737530
## iter   5 value -0.737530
## final  value -0.737530 
## converged

#ARMA_fit$ttable

Cuando analizamos graficamente los residuos vemos que en el modelo MA1 vemos: - Los residuos parecen tener un patròn - El ACF parece tener grandes valores en varios Lags - El Q-Statics muestra que todos los valores entan por debajo de la linea

Cuando analizamos gráficamente los residuos vemos que en el modelo MA2 vemos una mejora significtiva y posterioremente pasando al modelo ARMA vemos como mejoran aun mas, se reducen los valores del ACF, en el qqplot parecieran alinearse mejor los quantiles a la normalidad y asi tambien cuando analizamos el Q-statics los mismos obtienen un p-value mayor.

Analizaremos ahora el BIC y AIC de cada modelos, donde podemos observar que efectivamente y correspondiendose con lo observado hasta ahora el BIC y AIC mas bajo lo obtenemos con el modelo ARMA(1,1)

#Calculamos el AIC y BIC para cada uno de los modelos anteriores 
MA_fit$AIC
## [1] 1.401826
MA_fit$BIC
## [1] 1.422918
MA2_fit$AIC
## [1] 1.379757
MA2_fit$BIC
## [1] 1.40788
ARMA_fit$AIC
## [1] 1.375456
ARMA_fit$BIC
## [1] 1.403579

ARIMA

Hasta este momento trabajamos con modelos ARMA, introduciremos ahora el termino “I” por Integrated.

Para que ARIMA funcione al máximo, necesita que los datos estén estacionarios, esto es que la media y la varianza sean constantes en todo el conjunto. En los modelos se la denota con el termino “d”. Entonces, la diferenciación se utiliza para transformar los datos de modo que sean estacionarios.

Los modelos ARIMA como venimos viendo se definen por los ordenes (p,d,q).

Generamos a continuacion un modelo del orden 1,1,0, graficaremos posteriormente el ACF y PACF de los datos generados para ver cómo se comportan los datos integrados. Luego, diferenciaremos los datos para hacerlos estacionarios y comparemos su ACF y PACF.

x <- arima.sim(model = list(order = c(1, 1, 0), ar = .9), n = 200)
# Plot x
plot(x)

# Plot de ACF y PACF
acf2(x)

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  0.99  0.99  0.98  0.97  0.95  0.94  0.92  0.91  0.89  0.87  0.85  0.82
## PACF 0.99 -0.17 -0.13 -0.12 -0.10 -0.09 -0.07 -0.07 -0.06 -0.06 -0.05 -0.04
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.80  0.78  0.75  0.73  0.70  0.68  0.65  0.63   0.6  0.57  0.55  0.52
## PACF -0.02 -0.02  0.01  0.00 -0.01  0.00  0.00 -0.01   0.0 -0.01 -0.02 -0.03
##      [,25]
## ACF   0.50
## PACF -0.01
# Plot DIFF
plot(diff(x))

# Plot  PACF Y ACF 
acf2(diff(x))

##      [,1] [,2] [,3] [,4] [,5]  [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.92 0.85 0.78 0.73 0.69  0.64 0.62 0.61  0.57  0.53  0.51  0.48  0.44
## PACF 0.92 0.04 0.01 0.03 0.05 -0.01 0.08 0.10 -0.11 -0.05  0.07 -0.04 -0.09
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.40  0.36  0.33  0.28  0.24  0.19  0.18  0.16  0.13  0.10  0.07  0.03
## PACF  0.02 -0.03 -0.06 -0.12  0.03 -0.08  0.11  0.02 -0.10  0.02 -0.05 -0.09

Trabajaremos ahora con la base de datos de ASTSA globtemp que contiene desviaciones medias globales de la temperatura terrestre-oceánica, medidas en grados centígrados, para los años 1880-2015.

?globtemp
data("globtemp") #cargo la base de datos

Tanto el ACF como el PACF están disminuyendo, lo que implica un modelo ARIMA (1,1,1). El ACF se corta en el lag 2 y el PACF se está reduciendo, lo que implica un modelo ARIMA (0,1,2). El ACF se está reduciendo y el PACF se corta en el lag 3, lo que implica un modelo ARIMA (3,1,0). Aunque este modelo encaja razonablemente bien, puede ser el peor de los tres porque utiliza demasiados parámetros para autocorrelaciones tan pequeñas (como se puede ver cuando se analiza el BIC y AIC debajo)

# Plot  PACF y ACF de los datos 
acf2(diff(globtemp))

##       [,1]  [,2]  [,3] [,4]  [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.24 -0.19 -0.08 0.20 -0.15 -0.03  0.03 0.14 -0.16  0.11 -0.05  0.00
## PACF -0.24 -0.26 -0.23 0.06 -0.16 -0.09 -0.05 0.07 -0.09  0.11 -0.03 -0.02
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## ACF  -0.13  0.14 -0.01 -0.08     0  0.19 -0.07  0.02 -0.02  0.08
## PACF -0.10  0.02  0.00 -0.09     0  0.11  0.04  0.13  0.09  0.08
# Fit modelo ARIMA(1,1,1) 
ARIMA111_fit<-sarima(globtemp,1,1,1)
## initial  value -2.218917 
## iter   2 value -2.253118
## iter   3 value -2.263750
## iter   4 value -2.272144
## iter   5 value -2.282786
## iter   6 value -2.296777
## iter   7 value -2.297062
## iter   8 value -2.297253
## iter   9 value -2.297389
## iter  10 value -2.297405
## iter  11 value -2.297413
## iter  12 value -2.297413
## iter  13 value -2.297414
## iter  13 value -2.297414
## iter  13 value -2.297414
## final  value -2.297414 
## converged
## initial  value -2.305504 
## iter   2 value -2.305800
## iter   3 value -2.305821
## iter   4 value -2.306655
## iter   5 value -2.306875
## iter   6 value -2.306950
## iter   7 value -2.306955
## iter   8 value -2.306955
## iter   8 value -2.306955
## final  value -2.306955 
## converged

#ARIMA111_fit$ttable
# Fit aARIMA(0,1,2) 
ARIMA012_fit<-sarima(globtemp,0,1,2)
## initial  value -2.220513 
## iter   2 value -2.294887
## iter   3 value -2.307682
## iter   4 value -2.309170
## iter   5 value -2.310360
## iter   6 value -2.311251
## iter   7 value -2.311636
## iter   8 value -2.311648
## iter   9 value -2.311649
## iter   9 value -2.311649
## iter   9 value -2.311649
## final  value -2.311649 
## converged
## initial  value -2.310187 
## iter   2 value -2.310197
## iter   3 value -2.310199
## iter   4 value -2.310201
## iter   5 value -2.310202
## iter   5 value -2.310202
## iter   5 value -2.310202
## final  value -2.310202 
## converged

#ARIMA012_fit$ttable
# Fit ARIMA(3,1,0) 
ARIMA310_fit<-sarima(globtemp,3,1,0)
## initial  value -2.215090 
## iter   2 value -2.289035
## iter   3 value -2.306884
## iter   4 value -2.308838
## iter   5 value -2.309367
## iter   6 value -2.309746
## iter   7 value -2.309749
## iter   7 value -2.309749
## iter   7 value -2.309749
## final  value -2.309749 
## converged
## initial  value -2.314672 
## iter   2 value -2.314677
## iter   3 value -2.314679
## iter   4 value -2.314682
## iter   4 value -2.314682
## iter   4 value -2.314682
## final  value -2.314682 
## converged

#ARIMA012_fit$ttable
ARIMA012_fit$AIC
## [1] -1.723268
ARIMA012_fit$BIC
## [1] -1.637185
ARIMA111_fit$AIC
## [1] -1.716773
ARIMA111_fit$BIC
## [1] -1.630691
ARIMA310_fit$AIC
## [1] -1.717413
ARIMA310_fit$BIC
## [1] -1.60981

Forecasting ARIMA

Una vez que se elige un modelo, la proyecciòn es fácil.

Porque el modelo describe cómo se comporta la dinámica de la serie temporal a lo largo del tiempo. Simplemente continúa la dinámica del modelo en el futuro. En astsa, hay un comando llamado sarima-dot-for que se puede usar para pronosticar. Es similar al comando sarima, pero también especifica el horizonte de pronóstico.

Trabajemos con la base “oil” de la libreria ASTSA. Esta base de datos contiene datos de petróleo crudo, precio al contado del WTI FOB (en dólares por barril), datos semanales desde 2000 hasta mediados de 2010.

data("oil")
# ?oil

Usamos la función window para extraer subconjuntos de las series de tiempo. Para pronosticar los datos otro año, usamos sarima-dot-for y especificamos con qué anticipación deseamos pronosticar. Pronosticamos un avance de 52 semanas.

sarima.for imprime los pronósticos y sus errores estándar.

En el gráfico debajo las últimas 100 observaciones están trazadas en negro con puntos y los pronósticos están trazados en rojo. La muestra de color gris oscuro indica más o menos 1 error de predicción del cuadrado medio de la raíz. La muestra de color gris representa el intervalo de confianza del 95%. Además, el grafico incluye los precios reales del petróleo para 2007 (llamado oilf - f para el futuro) para que los valores predichos se puedan comparar con la verdad.

oil  <- window(astsa::oil, end = 2006) 
oilf <- window(astsa::oil, end = 2007)
sarima.for(oil, n.ahead = 52, 1, 1, 1)
## $pred
## Time Series:
## Start = c(2006, 2) 
## End = c(2007, 1) 
## Frequency = 52 
##  [1] 60.71882 60.43909 60.74214 60.75455 60.91190 60.99696 61.11808 61.22122
##  [9] 61.33332 61.44095 61.55081 61.65956 61.76887 61.87790 61.98706 62.09616
## [17] 62.20529 62.31441 62.42353 62.53265 62.64177 62.75089 62.86001 62.96913
## [25] 63.07825 63.18737 63.29649 63.40561 63.51473 63.62385 63.73297 63.84209
## [33] 63.95121 64.06033 64.16945 64.27857 64.38769 64.49681 64.60593 64.71505
## [41] 64.82417 64.93329 65.04241 65.15153 65.26065 65.36977 65.47889 65.58801
## [49] 65.69713 65.80625 65.91537 66.02449
## 
## $se
## Time Series:
## Start = c(2006, 2) 
## End = c(2007, 1) 
## Frequency = 52 
##  [1]  1.430557  2.270997  2.776644  3.245575  3.636008  3.996918  4.323903
##  [8]  4.629671  4.915600  5.186193  5.443159  5.688620  5.923876  6.150160
## [15]  6.368398  6.579407  6.783853  6.982317  7.175292  7.363212  7.546454
## [22]  7.725351  7.900198  8.071258  8.238767  8.402937  8.563961  8.722013
## [29]  8.877251  9.029820  9.179855  9.327476  9.472797  9.615922  9.756948
## [36]  9.895965 10.033055 10.168297 10.301764 10.433524 10.563640 10.692173
## [43] 10.819180 10.944712 11.068821 11.191554 11.312955 11.433067 11.551931
## [50] 11.669584 11.786062 11.901401
lines(oilf)

Hagamos ahora la proyecciòn de la temperatura global (globtemp ya cargada anteriormente) hasta el 2060, 45 años en el futuro. Recordemos que el modelo que hemos armado y que mejor se ajustaba a los datos era ARIMA(0,1,2)

ARIMA012_fit <- sarima(globtemp,0,1,2) # Fit ARIMA(0,1,2) 
## initial  value -2.220513 
## iter   2 value -2.294887
## iter   3 value -2.307682
## iter   4 value -2.309170
## iter   5 value -2.310360
## iter   6 value -2.311251
## iter   7 value -2.311636
## iter   8 value -2.311648
## iter   9 value -2.311649
## iter   9 value -2.311649
## iter   9 value -2.311649
## final  value -2.311649 
## converged
## initial  value -2.310187 
## iter   2 value -2.310197
## iter   3 value -2.310199
## iter   4 value -2.310201
## iter   5 value -2.310202
## iter   5 value -2.310202
## iter   5 value -2.310202
## final  value -2.310202 
## converged

ARIMA012_fit$ttable
##          Estimate     SE t.value p.value
## ma1       -0.3984 0.0808 -4.9313  0.0000
## ma2       -0.2173 0.0768 -2.8303  0.0054
## constant   0.0072 0.0033  2.1463  0.0337
#Hacemos la proyeccion 45 años para adelante desde el termino de la base de datos
sarima.for(globtemp, n.ahead=45, p=0, d=1, q=2) 
## $pred
## Time Series:
## Start = 2016 
## End = 2060 
## Frequency = 1 
##  [1] 0.7995567 0.7745381 0.7816919 0.7888457 0.7959996 0.8031534 0.8103072
##  [8] 0.8174611 0.8246149 0.8317688 0.8389226 0.8460764 0.8532303 0.8603841
## [15] 0.8675379 0.8746918 0.8818456 0.8889995 0.8961533 0.9033071 0.9104610
## [22] 0.9176148 0.9247687 0.9319225 0.9390763 0.9462302 0.9533840 0.9605378
## [29] 0.9676917 0.9748455 0.9819994 0.9891532 0.9963070 1.0034609 1.0106147
## [36] 1.0177685 1.0249224 1.0320762 1.0392301 1.0463839 1.0535377 1.0606916
## [43] 1.0678454 1.0749992 1.0821531
## 
## $se
## Time Series:
## Start = 2016 
## End = 2060 
## Frequency = 1 
##  [1] 0.09909556 0.11564576 0.12175580 0.12757353 0.13313729 0.13847769
##  [7] 0.14361964 0.14858376 0.15338730 0.15804492 0.16256915 0.16697084
## [13] 0.17125943 0.17544322 0.17952954 0.18352490 0.18743511 0.19126540
## [19] 0.19502047 0.19870459 0.20232164 0.20587515 0.20936836 0.21280424
## [25] 0.21618551 0.21951471 0.22279416 0.22602604 0.22921235 0.23235497
## [31] 0.23545565 0.23851603 0.24153763 0.24452190 0.24747019 0.25038375
## [37] 0.25326381 0.25611148 0.25892783 0.26171387 0.26447057 0.26719883
## [43] 0.26989951 0.27257344 0.27522139
lines(globtemp) #graficamos

Seasonal y Mix Models (SARIMA)

Trabajaremos ahora con modelos estacionales, si recordamos de cuando iniciamos estos videos y el inicio de esta notebook vimos por ejemplo el comportamiento de las ganancias trimestrales de Johnson & Johnson. Donde podiamos ver 1 ciclo por trimestre.

En los modelos SARIMA, la S es por seasonal.

Las series de tiempo puramente estacionales son raras y, por lo general, tenemos que mezclar la parte estacional con la parte no estacional.

Estos modelos mixtos se denominan SARIMA (p, d, q) x (P, D, Q) _S, las letras minúsculas denotan los órdenes de los componentes no estacionales y las letras mayúsculas se refieren a los componentes estacionales.

Trabajaremos sobre los datos de desemplo de EEUU por el periodo 1948-1978, de la libreria astsa que incluye los datos mensuales de desempleo de EE. UU.

Primero graficamos los datos para ver tendencias y persistencias estacionales. Luego observamos los datos sin tendencia, que deberian parecer estacionales.

data("unemp") #cargo monthly US unemployment data
#?unemp
# Plot unemp 
plot(unemp)

# Difference y grafico
d_unemp <- diff(unemp)
plot(d_unemp)

Ahora que hemos eliminado la tendencia y la variación estacional del desempleo, los datos parecen estacionarios.

# Seasonally difference d_unemp and plot it
dd_unemp <- diff(d_unemp, lag = 12)  
plot(dd_unemp)

Ajustamos un modelo SARIMA y observemos la muestra de ACF y PACF de la serie totalmente diferenciada.

Aclaracion: PACF y ACF esta en terminos anuales 1 año (12 meses), 2 años (24 meses), etc.

Vemos que en el componente no estacional el PACF se corta en el lag 2 y el ACF se reduce. Por otro lado en el componente estacional vemos que el ACF se corta en el lag 12 y el PACF se reduce en los lags 12, 24, 36,

# Plot PACF y ACF con lag 60
dd_unemp <- diff(diff(unemp), lag = 12)
acf2(dd_unemp, max.lag = 60)

##      [,1] [,2] [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.21 0.33 0.15 0.17 0.10  0.06 -0.06 -0.02 -0.09 -0.17 -0.08 -0.48 -0.18
## PACF 0.21 0.29 0.05 0.05 0.01 -0.02 -0.12 -0.03 -0.05 -0.15  0.02 -0.43 -0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.16 -0.11 -0.15 -0.09 -0.09  0.03 -0.01  0.02 -0.02  0.01 -0.02  0.09
## PACF  0.15  0.03 -0.04 -0.01  0.00  0.01  0.01 -0.01 -0.16  0.01 -0.27  0.05
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF  -0.05 -0.01  0.03  0.08  0.01  0.03 -0.05  0.01  0.02 -0.06 -0.02 -0.12
## PACF -0.01 -0.05  0.05  0.09 -0.04  0.02 -0.07 -0.01 -0.08 -0.08 -0.23 -0.08
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
## ACF   0.01 -0.03 -0.03 -0.10 -0.02 -0.13  0.00 -0.06  0.01  0.02  0.11  0.13
## PACF  0.06 -0.07 -0.01  0.03 -0.03 -0.11 -0.04  0.01  0.00 -0.03 -0.04  0.02
##      [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
## ACF   0.10  0.07  0.10  0.12  0.06  0.14  0.05  0.04  0.04  0.07 -0.03
## PACF  0.03 -0.05  0.02  0.02 -0.08  0.00 -0.03 -0.07  0.05  0.04 -0.04
# Fit del modelo
sarima(unemp, p = 2, d = 1, q = 0, P = 0, D = 1, Q = 1, S = 12)
## initial  value 3.340809 
## iter   2 value 3.105512
## iter   3 value 3.086631
## iter   4 value 3.079778
## iter   5 value 3.069447
## iter   6 value 3.067659
## iter   7 value 3.067426
## iter   8 value 3.067418
## iter   8 value 3.067418
## final  value 3.067418 
## converged
## initial  value 3.065481 
## iter   2 value 3.065478
## iter   3 value 3.065477
## iter   3 value 3.065477
## iter   3 value 3.065477
## final  value 3.065477 
## converged

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     ar2     sma1
##       0.1351  0.2464  -0.6953
## s.e.  0.0513  0.0515   0.0381
## 
## sigma^2 estimated as 449.6:  log likelihood = -1609.91,  aic = 3227.81
## 
## $degrees_of_freedom
## [1] 356
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1    0.1351 0.0513   2.6326  0.0088
## ar2    0.2464 0.0515   4.7795  0.0000
## sma1  -0.6953 0.0381 -18.2362  0.0000
## 
## $AIC
## [1] 8.991114
## 
## $AICc
## [1] 8.991303
## 
## $BIC
## [1] 9.034383

Forecasting SARIMA

Proyectemos ahora el desemplo para los proximos 3 años (36 meses).

sarima.for(unemp, n.ahead=36, p=2,d=1, q=0, P=0, D=1,Q=1,S=12)
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1979 676.4664 685.1172 653.2388 585.6939 553.8813 664.4072 647.0657 611.0828
## 1980 683.3045 687.7649 654.8658 586.1507 553.9285 664.1108 646.6220 610.5345
## 1981 682.6406 687.0977 654.1968 585.4806 553.2579 663.4398 645.9508 609.8632
##           Sep      Oct      Nov      Dec
## 1979 594.6414 569.3997 587.5801 581.1833
## 1980 594.0427 568.7684 586.9320 580.5249
## 1981 593.3713 568.0970 586.2606 579.8535
## 
## $se
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1979  21.20465  32.07710  43.70167  53.66329  62.85364  71.12881  78.73590
## 1980 116.99599 124.17344 131.51281 138.60466 145.49706 152.12863 158.52302
## 1981 194.25167 201.10648 208.17066 215.11503 221.96039 228.64285 235.16874
##            Aug       Sep       Oct       Nov       Dec
## 1979  85.75096  92.28663  98.41329 104.19488 109.67935
## 1980 164.68623 170.63839 176.39520 181.97333 187.38718
## 1981 241.53258 247.74268 253.80549 259.72970 265.52323
lines(unemp) #graficamos

Bibliografia

El contenido de esta notebook resulta de una mezcla de distintos lugares, entre otros: